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We obtain a convergent power series expansion for the first branch of the dispersion 
relation for subwavelength plasmonic crystals consisting of plasmonic rods with 
frequency-dependent dielectric permittivity embedded in a host medium with unit 
permittivity. The expansion parameter is n = kd — 2ird/\, where k is the norm 
of a fixed wavevector, d is the period of the crystal and A is the wavelength, and 
the plasma frequency scales inversely to d, making the dielectric permittivity in the 
rods large and negative. The expressions for the series coefficients (a.k.a., dynamic 
correctors) and the radius of convergence in r\ arc explicitly related to the solutions 
of higher-order cell problems and the geometry of the rods. Within the radius of 
convergence, we are able to compute the dispersion relation and the fields and 
define dynamic effective properties in a mathematically rigorous manner. Explicit 
error estimates show that a good approximation to the true dispersion relation is 
obtained using only a few terms of the expansion. The convergence proof requires 
the use of properties of the Catalan numbers to show that the series coefficients are 
exponentially bounded in the H 1 Sobolcv norm. 

Keywords: Meta-material, Plasmonic crystal, Dispersion relation, Effective 
property, Series solution, Catalan number 

1. Introduction 

Sub-wavelength plasmonic crystals are a class of meta-material that possesses a 
microstructure consisting of a periodic array of plasmonic inclusions embedded 
within a dielectric host. The term "sub- wavelength" refers to the regime in which 
the period of the crystal is smaller than the wavelength of the electromagnetic 
radiation traveling inside the crystal. Many recent investigations into the behav- 
ior of meta-materials focus on phenomena associated with the quasi-static limit in 
which the ratio of the period cell size to wavelength tends to zero. Sub-wavelength 
micro-structured composites are known to exhibit effective electromagnetic proper- 
ties that are not available in naturally-occurring materials. Investigations over the 
past decade have explored a variety of meta-materials, including arrays of micro- 
resonators, wires, high-contrast dielectrics, and plasmonic components. The first 
two, especially in combination, have been shown to give rise to unconventional bulk 
electromagnetic response at microwave frequencies (Smith et al. 2000; Pendry et 
al. 1999; Pendry et al. 1998) and, more recently, at optical frequencies Povinclli 
(2009), including negative effective dielectric permittivity and/or negative effective 
magnetic permittivity. An essential ingredient in creating this response are local 
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0.1 


0.2 


0.3 


0.4 


0.45 


Rm 


1/60 


1/68 


1/88 


1/96 


1/340 



Table 1. Lower bounds on the radii of convergence R for circular inclusions of radii rd. 
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0.1 


0.2 


0.3 


0.4 


0.45 


A m 


38 pan 


43 fim 


56 fim 


73 fim 


214/xra 




1.6 • 10 5 m _1 


1.4 • 10 5 m _1 


i.i ■ w 6 m - 1 


1.0- 10 5 m" 1 


3.0 • 10 4 m" 1 



Table 2. Values of A m and km for circular inclusions of radii rd when d = 10 7 m. 



resonances contained within each period due to extreme properties, such as high 
conductivity and capacitance in split-ring resonators Pendry et al. (1999). 

In the case of plasmonic crystals, the dielectric permittivity e p of the inclusions 
is frequency dependent and negative for frequencies below the plasma frequency lo p , 

e P H = l-4- (1-1) 

Shvets & Urzhumov (2004, 2005) have investigated plasmonic crystals in which uj p 
is inversely proportional to the period of the crystal and for which both inclusion 
and host materials have unit magnetic permeability. They have proposed that si- 
multaneous negative values for both an effective e and [i arise at sub-wavelength 
frequencies that are quite far from the quasi-static limit, that is, 

2?rd , N 

r] = kd=— (1.2) 
A 

is not very small, where d is the period of the crystal, k is the norm of the Bloch 
wavevector and A is the wavelength. In this work, we present rigorous analysis of this 
type of plasmonic crystal by establishing the existence of convergent power series 
in r\ for the electromagnetic fields and the first branch of the associated dispersion 
relation. The effective permittivity and permeability defined according to Pendry 
et al. (1999) are shown to be positive for all rj within the radius of convergence 
R, and, in this regime, the extreme property of the plasma produces no resonance 
in the effective permittivity or permeability. This regime is well distanced from the 
resonant regime investigated in Shvets & Urzhumov (2004, 2005). 

The analysis shows that the radii of convergence of the power series is at least 
R m , which is not too small, as shown in Table [l] which contains values of R m for 
circular inclusions of various radii rd. The number R m can be put in physical per- 
spective by fixing the cell size and introducing the parameters A m and &m such that 
the power series describes wave propagation for wavelengths above X m and wave 
numbers below km- Table 2 contains values of A m and &m when d = 10~ 7 to. The 
wavelengths lie in the infrared range and the plasma frequency is uj p — 10 15 sec _1 . 

We focus on harmonic H-polarized electromagnetic waves in a lossless compos- 
ite medium consisting of a periodic array of plasmonic rods embedded in a non- 
magnetic frequency-independent dielectric host material. Each period can contain 
multiple parallel rods with different cross-sectional shapes, however the rod-host 
configurations are restricted to those with rectangular symmetry, i.e., configura- 
tions invariant under a 180° rotation about the center of the unit cell. The regime 
of interest for this investigation is that in which 
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Figure 1. Unit cell with plasmonic inclusion. 



1. the plasma frequency lo v is high 

2. the ratio of the cell width to the wavelength is small (rj <C 1). 

From the formula e p = 1 — ui 2 /ui 2 , it is seen that a high plasma frequency ui p gives 
rise to a large and negative dielectric permittivity e p in the plasmonic inclusions. 
Following Shvets & Urzhumov (2004), the plasma frequency is related to the cell 
size by 



This results in the relation 



1 



where c is the speed of light in vacuum. The governing family of differential equa- 
tions for the magnetic field is the Helmholtz equation with a rapidly oscillating 
coefficient 

LO 2 

-V ■(A(x/d)Vu) = -5-u, (1.3) 
& 

in which A is the matrix defined on the unit period of the crystal by 

. . . f e p -1 / in the plasmonic phase, 
1 c^ 1 ! in the host phase, 

/ is the identity matrix and 



1 and e„ = 1 



? -d 2 ' 



The coefficient A is not coercive in the regime lu v > u>, since e p is negative in this 
regime, and it is precisely the appearence of negative e p that allows us to obtain a 
convergent power series expansion of the electromagnetic field and the frequency u> 2 
for a fixed Bloch wavevector k = kk, with \k\ = 1. 



In Theorem (5.2 1, we obtain the following series expansion for the frequency u> 2 

oo 

u 2 = u 2 = c 2 k 2 £ £ m V 2m (1.4) 

m=0 

in which £ 2 m is a tensor of degree 2m + 2 in k. This gives rise to a convergent power 
series for an effective index of refraction n 2 s defined through 

r 2 k 2 

n eff=— (I- 5 ) 
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The effective property n^ s is well defined for all r\ in the radius of convergence and 
is not phenomenological in origin but instead follows from first principles using the 
power series expansion. Interpreting the first term of this series as the quasi-static 
index of refraction n^, the remaining terms then provide the dynamic correctors 
of all orders. In section 6, we define the effective permeability fi c s and prove that 
n^ ff and /j, e g are both positive for r\ in some interval (0, tjq] and that a mild effective 
magnetic response emerges for the homogenized composite, even though the com- 
ponent materials are non-magnetic (/ip = /ip = 1). Having defined n^ ff and /x e ff, 
the effective electrical permittivity e c g can be defined through the equation 

2 

so that e c ff is positive whenever both and /i c ff are positive. Thus, one has 
a solid basis on which to assert that plasmonic crystals function as materials of 
positive index of refraction in which both the effective permittivity and perme- 
ability are positive. The method developed here can be applied to other types of 
frequency-dependent dielectric media such as polaratonic crystals. From a physical 
perspective, this work provides the first explicit description of Bloch wave solutions 
associated with the first propagation band inside nanoscale plasmonic crystals. In 
the context of frequency independent dielectric inclusions, the first two terms of 
n^ ff are identified via Rayleigh sums in McPhedran et al. (2006). 

To emphasize the difference between effective properties defined for meta-material 
structures where the crystal period d is fixed and effective properties defined in the 
quasistatic limit, i.e., k fixed and d — > 0, we refer to the latter as quasistatic effective 
properties and denote these with the subscript qs. The situation considered in this 
paper contrasts with the case in which e ~ d~ 2 in the inclusion and is large and 
positive, investigated by Bouchitte & Felbacq (2004). In that case for 77 — i> 0, n qs (io) 
has poles at Dirichlct eigenvalues of the inclusion and therefore is negative in cer- 
tain frequency intervals (see also Bouchitte & Felbacq (2004), (2005), (2005a)). In 
fact, what allows us to prove convergence of the power series in the plasmonic case 
is precisely the absence, due to negative e p , of these internal Dirichlet resonances. 

In the regime where e p is negative and large, the perturbation methods used for 
describing Bloch waves in heterogeneous media developed in Odeh & Keller (1964), 
Conca (2006), Bensoussan et al. (1978) cannot be applied. Our analysis instead 
makes use of the fact that e p is negative and large for sub-wavelength crystals and 
develops high-contrast power series solutions for the nonlinear eigenvalue problem 
that describes the propagation of Bloch waves in plasmonic crystals. The conver- 
gence analysis takes advantage of the iterative structure appearing in the series 
expansion and is inspired by a technique of Bruno (1991) developed for series solu- 
tions to quasi-static field problems. We prove that the series converges to a solution 
of the harmonic Maxwell system for ratios of cell size to wavelength that are not 
too small. Indeed for typical values of the plasma frequency w p the analysis delivers 
convergent series solutions for nano scale plasmonic rods at infrared wavelengths. 

In section [6] we compute the first two terms of the dispersion relation for circular 
inclusions Shvets & Urzhumov (2004, 2005) and provide explicit bounds on the 
relative error comitted upon replacing the full series with its first two terms. The 
error is seen to be less than 3% for values of rj up to 20% of the convergence 
radius, so that the two-term approximation provides a numerically fast and accurate 
approximation to the dispersion relation. 
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The high contrast in e gives rise to effective constants e e g and /j, e g. In the bulk 
relation 

B c g = A'cff-Hcff, (1-6) 

where B c s is the average over the the period cell (a flux), whereas H e g is the 
average of H3 over line segments in the matrix parallel to the rods. Taking the 
ratio of B c ff/H c ^ delivers an effective magnetic permeability and one recovers mag- 
netic activity from meta-materials made from non-magnetic materials. This phe- 
nomenon was understood by Pendry et al. (f999), and has been made rigorous in 
the quasistatic limit through two-scale analysis in several cases. These include the 
two-dimensional arrays of inclusions in which e scales as dr 2 (Bouchitte & Felbacq 
(2004, 2005); Felbacq & Bouchitte (2005, 2005a); Pendry et al. (1999)) two dimen- 
sional arrays of ring resonators whose surface conductivity scales as d _1 Kohn & 
Shipman (2008), as well as three-dimensional arrays of split-ring wire resonators 
in which the conductivity scales as dr 2 Bouchitte & Schweizer (2008). This "non- 
standard" homogenization has been understood for some decades in problems of 
porous media and imperfect interface (Cioranescu & Paulin (1979); Auriault & Ene 
(1994); Lipton (1998); Donnato & Monsouro (2004); Zhikov (1995)) and recently 
has given rise to interesting effects in composites of both high contrast and high 
anisotropy (Cherednichcnko et al. (2006); Smyshlyaev (2009)). 

The two-scale analysis in these cases relies on the coercivity of the underlying 
partial differential equations. The problem of plasmonic inclusions, however, is not 
coercive because e is negative in the plasma — but it is precisely this negative index 
that underlies the convergence of the power series. As we shall see, the uniqueness of 
the solution of the Dirichlet problem for Au — u = in the plasmonic inclusion gives 
exponential bounds on the coefficients of the series, which allows us to prove that 



it converges to a solution of the differential equation (1.3). This result presented 
here shows that by considering a finite number of terms in the series, one has an 
approximation of the true solution, to any desired algebraic order of convergence. 
In this context we point out the recent work of (Smyshlyaev & Cherednichcnko 
(2000); Kamotski et al. (2007); Panasenko (2009)) that shows that the power series 
expressed by the formal two scale expansion of Bakhvalov & Panasenko (1989) is 
an asymptotic series in certain cases under the hypothesis that the coefficient A is 
coercive. 



2. Mathematical Formulation and Background 

We introduce the nonlinear eigenvalue problem describing the propagation of Bloch 
waves inside a plasmonic crystal and provide the context for the power series ap- 
proach to its solution. 

For points x = (21,2:2) in the sia^-plane, the d-periodic dielectric coefficient of 
the crystal is denoted by e(w,x), where 



e(w,x) 



e p (w) for x e P, 
e 5 for x G P. 



Both materials are assumed to have unit magnetic permeability, /i p = /ip = 1. 

We assume a Bloch- wave form of the field, where k = (ki, K2) is the un it vector 
along the direction of the traveling wave and k — 2ir/\ is the wave number for a 
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wave of length A. The magnetic and electric fields are denoted by H = (Hi, H2, H 3 ) 
and E = (E±, E 2 , E 3 ) respectively. For P-polarized time-harmonic waves, the non- 
vanishing field components are 

H 3 = H a (x)e ilk *' x - U *\ E x = Pi(x)e^ x - Wt) , E 2 = E 2 {x)e r{kk ^-^ (2.1) 

in which the fields H 3 (x), Pi(x), and ^(x) are continuous and d-periodic in both 



X\ and X2- The Maxwell equations take the form of the Helmholtz equation (1.3 1 
in which substitution of u = H 3 (x.)e ll - kK ' x ^ ut ' > gives 



uJ 2 



(V + ikk)e~ L (V + ikk)H 3 = -^H 3 in the rods, (2.2) 

c 2 



(V + ikk)e- 1 (V + ikk)H 3 — in the host material, (2.3) 



c 2 



where H 3 satisfies the transmission conditions on the interface between the rods 
and host material given by 

n-(e p -\V + ikk)H 3 ) lp =n-(e^(V + ikk)H 3 ) lp . (2.4) 

Here, the subscripts indicate the side of the interface where the quantities are 
evaluated and n is the unit normal vector to the interface pointing into the host 
material. We denote the unit vector pointing along the x 3 direction by e3, and the 
electric field component of the wave is given by E = — ^e 3 x VH 3 . 



For each value of the wave- vector k, equations (2.2 2.3 2.4 1 provide a nonlinear 
eigenvalue problem for the pair H 3 and lo 2 . One of the main results of this work is 
to show that this problem is well posed by explicitly constructing solutions using 
power series expansions. In order to develop the appropriate expansions, we rewrite 



(2.2 2.3 2.4 1 in terms of 77 and a dimensionless variable y in M 2 that normalizes a 
period cell to the unit square Q = [0, l] 2 , x = yd = yrj/k. We define the Q-periodic 
function 

My) = H 3 (yd) 
and for convenience of notation, we redefine 

e ( W ,y) = j £pM f ryG |' (2.5) 
V 1 e- p foryeP, y ' 

to arrive at the eigenvalue problem that requires the pair h(y) and w 2 to be a 
solution of the master system 

-(V y + «^)r 1 (w,y)(V y + j»jK)/i(y) = ?7 2 ^My) for ye PUP, 
n ■ e p " 1 (w)(V y + ir/k)h(y)\ p = n ■ ef^Vy + ir]k)h(y)\p for y G dP. 

We prove in Theorem |5.2| that this eigenvalue problem can be solved by constructing 
explicit convergent power series solutions. 

The development of the remainder of the paper is as follows. In section [3] the 
power series expansion is introduced and the associated boundary-value problems 
necessary for determining each term in the series are obtained. The boundary value 
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problems are given by a strongly coupled infinite system of linear partial differential 
equations. The existence and uniqueness of the solution to this infinite system is 
proved under fairly general hypotheses in section [4] Because the system is coupled 
through convolution products, the convergence analysis is delicate. The convolu- 
tions are handled through estimates involving sequences of Catalan numbers whose 
convolution products determine the next element of the sequence. The Catalan 
numbers and their relevant properties are discussed and used to derive bounds on 
the Sobolcv norm of each term of the series expansion in section [5j These bounds 
are then used to establish the radius of convergence for the power series representa- 
tions of the field and frequency (the dispersion relation) , which solve the nonlinear 
eigenvalue problem (2.6 1. Section [6] deals with the computation of error bounds for 



2 



finite-term approximations of the magnetic field and the dispersion relation. 

3. Power Series Expansions 

We take r\ to be the expansion parameter for the field h(y) and the frequency ui 

h v = h + 7]hi + rfh 2 + 

<4 = Uo+Vul + V 2 u% + ---, (3-1) 
in which the functions h m are periodic with period cell Q. 



Inserting (3.1 1 into (2.6 1 and identifying coefficients of like powers of 77 on the 



right- and left-hand sides yields the equations 

2 

Ah m + 2ik-Vh m ^ 1 - /i m „ 2 = -^2^m-2-£ in P, 

2 

A/i m + 2iK-V/i m _i - /i m _ 2 = h m - ^/i m - 2 -£ in P, 

(^prVh, m _2-^ - V/i m - ih m ^ 1 kj \ P ■ n = ^feVh m - 2 -t\ p ■ n on dP, 

for to = 0, 1,2,..., in which h m = and ujf n — for m < and the terms involv- 
ing the subscript t are convolutions written according to the following summation 
conventions, 



a e b n -i = ^ a tK-£ , a e b n-t 2) = T"! 

e=o 1=0 

e 2 -l [n/2] 

Ati<t<l2) ST^ u At even) ST^ 

i=£i+i e=o 



-21, 



where [n/2] denotes the largest integer less than or equal to n/2. The boundary 
value problem satisfied by Hq in P is 

Ah = in P, 

V/i |p-n = ondP, 

so that this function is necessarily a constant in P. We denote this constant value 
by Jiq. It will be convenient to use the dimensionless parameters ip m and defined 
through 

h m = i m h ip m and u 2 m = c 2 k 2 f m . (3.2) 
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In terms of ip m and the above equations for the functions h m become 

Alp m + 2k ■Vl/'m-l + 1pm-2 = (-iY^1p m -2-£ ™ P, 

Alp m + 2k -Vlprn-1 + lprn-2 = Ipm + {-iYQ^m-2-t ™ P, 

(V((-j)^^Vm-2-^) + V^m + ^m-l«)|j> ' n = V((-i)*£f )|p ' n On 9P, 

(3.3) 

for to = 0,1,2,..., in which = and £^ = for m < 0. We thus have 
an infinite system which the sequences {ip m }m=o ano - {£m}m=o must satisfy. The 
system is written in terms of a Poisson equation in P with Neumann boundary 
data and a Helmholtz equation in P with Dirichlet boundary data, namely 



and 



in which 



Ai/j m = G m in P, 

Vip m \p n = F m ondP, 

A^ m = ip m + G m in P, 

ip m \ P = ipm\p on dP, 



(3.4) 



(3.5) 



F m = V((-lYe^rn-2-e)\ P ■n-(V((-iye e ^rn-2-i) ~ 



(3.6) 



G m = {-iY£,}i> m -2-t - 2ik-\7ip m -! - tpm-2 
where ip m \p and ip m \ v indicate the trace of tp m on the P and P sides of the interface 



dP separating the two materials. The Neumann boundary value problem (3.4 1 is 
subject to the standard solvability condition given by 

(G m )p+(F m ) dP = 0. (3.7) 

Here the area integral over the domains P and P are denoted by (-)p and (-)p, while 
the line integral over the interface dP separating the two materials is denoted by 
(■)op. 

The iterative algorithm for solving the system is as follows. First note from the 
definition of ipo it follows that that ipo = I for y in P. The function i/jq is determined 



inside P by solving (3.5) with Dirichlet boundary data ^o| p = ipo\p = 1 on dP. 



Then ipi on P is the solution of (3.4) with Neumann data Vtpi\p ■ n = —ipo\pk ■ n 
on dP. The process then continues with the boundary values on dP of ip m in P 
providing the Dirichlet data for ip m in P which, in turn, provides the Neumann data 
for ipm+i in P, up to an additive constant. The term is determined by the 



consistency condition (3.7 1 and an inductive argument can be used to show that it 
is a monomial of degree m in k. The equations satisfied by -0 n , . . . , i/'4 inside P, P, 
and dP are listed in Table [3] below. 

Note in the table that £% dd — (meaning £| = for £ = 1, 3, 5, . . . ). In the next 
section we identify a large class of shapes for the plasmonic rod cross sec tions for 
which the sequences {ipm(y)}m=o ail d {£m}m=o satisfy the infinite system (3.4 3.5 



3.6 3.7 1 and (ipm) p = 0, to = 1, 2, . . .. The mean zero property of ip m on P provides 
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p 


P 


9F 


0o = i 


Ai/>o = i/'o 


VVolp - n = 


Ai/>i + 2k Vi/>0 = 


Ai/ji + 2k-V0o = -01 


(Wi + 0ok)Ip • n = 


Aj/> 2 +2ft-Wi + 0o = ?o0>o 


A^ 2 + 2ftVV>i + 0o = V>2 + SoV'o 


(V(^Vo) + V02 + 0i«)|p-n 
= V(^o)| P -n 


A7/J3 +2K-V02 +Vi = ?o0>i 


A0 3 + 2K-V02 + 01 = V>3 + 5oVl 


(V(^Vl)+V0 3 +02ft)|p-n 
= V(f^l)l P -n 


A</>4 + 2/t-Vl/> 3 + 02 = £o02 - glpO 


A0 4 + 2K-V03 + 02 
= (?O02~ gi>o) +04 


(V(^0 2 - + + V-3«)Ip • n 
= V(^V>2-«i0o)l P -n 



Table 3. Table of PDEs for the functions tp m obtained from the expansion in n. 



a tractable scenario for proving the convergence of the resulting power series. This 
topic is discussed further in section [4] 

In what follows we will make use of the equivalent weak form of the infinite 
system. To introduce the weak form we introduce the space of complex valued 
square integrable functions with square integrable derivatives H 1 (Q). For u and v in 

H l {Q) the inner product is defined by (u, w)_ff 1 (Q) = ( Sq uv dy + fg Vu • Vvdy 

and the norm is given by ||i>||.h-i(q) = ( v i v ) X h^(q)- ^ nc H 1 inner products and norms 
over P and P are defined similarly. 

The weak form of the infinite system is given in terms of the space Hp C1 .(Q) of 
functions in H 1 (Q) that take the same boundary values on opposite faces of Q. The 
weak form of the system (3.4 3.5 3.6 3.7 1 is given by 



+([Vf/>m + kipm-i] ■ Vw - [k -Vlpm-l + 1p m -2]v) p = Q, 



(3.8) 



for all v G Hl er (Q), where a' m = (-i)'£?Vm-i and o"m = (-i) e i>m-t£t 



,£ 2 



The 



equivalence between (3.4 3.5 3.6[ ) and the weak form follows from integration by 
parts and the solvability condition (13 .Th follows from (3.8 1 on choosing the test 



function v = 1 in (3.8 1. 



4. Solutions of the Infinite System for Plasmonic Domains 
with Rectangular Symmetry 

The goal here is to identify solutions of the infinite system for which one can prove 
convergence of the associated power series with a minimum of effort. Looking ahead 
we note that the convergence proof is expedited when one can apply the Poincare 
inequality to the restriction of tp m on P for m greater than some fixed value. To 
this end we seek a solution ({tp m (y)}m=o> {£m}m=o) sucn that for m > 1 one has 
(ipm)p — and the sequences satisfy satisfy ( 3.4| 3.5 3.6 3.7 1 or equivalently satisfy 



( 3.8 1 . We show that we can find such solutions for the class of plasmonic domains P 
with rectangular symmetry. Here we suppose that the unit period cell is centered at 
the origin and the class of rectangular symmetric domains is given by the set of all 
shapes invariant under 180° rotations about the origin. This class includes simply 
connected domains such as rectangles and ellipses as well as multiply connected 
domains. For these geometries and for each m = 1, 2, 3 ... it is demonstrated that 
one can add an arbitrary constant to the restriction of the function tp m on P with 



out affecting the solvability condition (3.7). Under the assumption of rectangular 
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symmetry for the inclusion P, we will show that there exists a pair ({VVn(y)}m=o> 
{£m}m=o) satisfying (3.8) with the functions ipm in the subspace H* (Q) C H* CY .(Q) 
of real-valued junctions with zero average in P. 

We now record the symmetries necessarily satisfied by any solution ip m £ Hl(Q) 



to (3.4 3.5 3.6 1 for plasmonic domains with rectangular symmetry. We denote the 
on the unit vector k by writing ip 1 ^ so that 



dependence of tp, 
(i) ^-«(y) = ( 



-irvd(2/),VyeQ 



(ii) ip-*(y) = iPU-y), Vy e Q. 

Statement (i) is true for inclusions of arbitrary shape, while statement (ii) is true 
only for inclusions with rectangular symmetry. Taken together, these statements 
imply that ip^—y) = (—l) m ^p^ l (y), so that i^m is even or °dd in Q according as 
the index m is even or odd. From its definition, tpQ = 1 in P and trivially satisfies 
the solvability condition (3.7). The solvability of ip m when m > 1 is proved by 



induction on to using the weak form (3.8 1. We have the following theorem 



Theorem 4.1. For each k, there exists a sequence of junctions {ip m }m=i> VVn € 
H}(Q), and a sequence of real numbers with £,g dd = 0, solving the weak form 

for each integer to. 

Proof. The proof is divided into the base case (to = 1 and m — 2) and the inductive 
step. 

Base case : 

The solvability for ipi and ip2 can be established without the need to restrict 
to rectangular symmetric inclusions. This restriction will be necessary only in the 
inductive step. Setting to — 1 and v = 1 in (3.8 1, we see that the left-hand side of 
(3.8 1 vanishes. This establishes the solvability for ip\. If we then take (ipi)p = 0, 
we have a solution ipi £ H}(Q). Setting to = 2 and v = 1 in (3.8), we obtain 



(a' ) P + (a' ) p - (k • V^i + Vo) P = 0. 

Since (V'o)q > (see Appendix) and (cr' }p + (cr' }p — £,o(^o)q, this is one equation 
in one unknown £q. Solving for Q then gives £q = (V'o)q 1 ('« • V^i + V-'o)p- Choosing 
this value for and also taking (ip2)p = 0, we have a solution ip 2 & Hl(Q). 
Inductive step: 

Let 2n be an even positive integer and assume that (3.8 1 has solutions ip m £ 
Hl(Q) for m = 1,2, ...,2n, with £_ 2 e R and £ 

-02«+i, '02n+2 € for m = 2n - 



</</ 



l,2n+ 2 with £| T 



0. Then (3.8 1 has solutions 



and ei» S 



The solvability condition for ip2n+i is obtained by setting w = 1 and m — 2n + 1 
in the weak form, namely 

([«' Vcr 2n-2 - °2n-l ~ °2n-3 + ^Lsl) P + 
+ ([« ; -Vo-2 n -2 - °2n-l - CT 2«-3 + ^2n-s])p + 
+{[K-V'02n + V ) 2n-l])p = O. 

The hypothesis = 0, odd < 2n — 2, will imply that the convolutions cr m , 
to < 2n — 2, have the same even/odd property as the functions ip m . Indeed, writing 
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out o"2 n _3, we have cr 2n „ 3 — (— i) e £,iip 2n ^_^, anc * since 2n — 3 — £ is odd when £ is 
even, it follows that <7 2 n-3 is a linear combination of odd functions and is, therefore, 
an odd function. The same reasoning applies to all the other convolutions of index 
less than or equal to 2n — 2. Moreover, K-\7a' 2n _ 2 is an °dd function, since cr' 2n _ 2 is 
even. Thus, all integrals in the consistency condition above vanish (for the integrals 
in P we can also use the fact that all functions belong to H}(Q)), except that 

«-i)p + «-i)p = (-^-^L-ii^Q. 

Since (V'o)q > (see Appendix), the solvability condition for V^n+i is simply 
Cln-i = 0- We thus take £fn-i = to establish the existence of ip2n+i = 0. 
Moreover, since ipm and £ m _2 are real by the induction hypothesis, < m < 2n, 
it follows that ip2n+i is real-valued. Thus, taking (tp 2n+ i)p — 0, we have a solution 
V-'2n+i € Hl(Q). Also, ip2n+i is an odd function since its index is odd. We now 
proceed to the solvability of ip2n+2 , namely 

([£-Vcr 2 n-l - 02n-2 - <4n-2 + <4J)p + 
+ <[K-Vcr 2 „_i - o' 2n _ 2 - 0- 2 '„_ 2 + CT 2n ])p + 
+ ([k -V-02«+l + ^2n]) P = 0. 

All terms in the above equation are real numbers, since we assumed ip m and 
£m-2 rea l f° r < m < 2n, with <^ dd = 0, and we just took £fn-l = ^ anc ^ 
"02n+i is real-valued. Thus, this equation contains the only one undetermined term 
(— i) 2n Sani^o) Q ■ Thus, we have one real equation with one real variable, so that 
taking £| Tl to be such as to solve this equation and also taking (ip2n+2)p = 0, we 
complete the proof of the inductive step. □ 



5. Proof of Convergence 

In this section we show that the power series J2m=oP mT l m > J2m=oPmV m and 
ELoCf. where p m = \\^ m \\m(P) and p m = \\ip m \\ H ^p), converge and pro- 
vide lower bounds on their radius of convergence. This will then be used to show 
that the pair h v = J2m=ohmV m and ui^ — ^m^o^mV" 1 is a solution to (2.6 1. In 



subsection [a] we present the Catalan Bound, which is used to provide a lower bound 
on the radius of convergence of the power series. In subsection [b] we derive inequal- 
ities which bound p m , p m and £^ in terms of lower index terms. In subsection [c] 
we present the properties of the Catalan numbers relevant for bounding convolu- 
tions of the kind appearing in (3.4 1 and (3.51. In subsection [d] we use an inductive 



argument on the inequalities of subsection [b] to prove the Catalan Bound. Finally, 
in subsection |e] we prove that the pair /i„ and ^ is a solution to the eigenvalue 
problem Kffy. 



(a) The Catalan Bound 

The following theorem is one of the central results of this paper 

Theorem 5.1. (Catalan Bound) 
For every integer m, we have that 

Pm, Pm, \C\ < PC m J m (5.1) 
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in which C m is the m th Catalan number, (3 = max{po, po, |£ 2 |} an d J — Tnax{J%, J2}, 
where the numbers J\ and Ji are determined as follows: J\ is the smallest value of 



J such that (5.1l holds for m < 4 and J2 is the smallest value of J for which the 
following polynomials Q* , R* , S* in the variable J^ 1 are all less than unity 

Q* = ftp[ A{2£;(4)/3J~ 2 l/3+£;(4)/3J" 3 5/42+£;(4)/3J~ 4 l/21+£; 2 (4)/3 2 J~ 4 5/42}- 
+ 2E(4)/3 J- 2 l/3 + 2E (A) (3 J- 3 5/ 42 + £(4)/3 J~ 4 1/21 + E 2 {4)0 2 J~ A b/42 + J~ 2 5/42 
+ 2n P (A{2E(4)f3J- 3 5/42 + 2E(4)l3J- 4 l/21 + E(4)l3J- 5 l/42+E 2 (4)f3 2 J- 5 l/21}+ 
+ 2£(4)/3J~ 3 5/42 + 2E(4)f3J- 4 l/21 + E(4)/3J~ 5 l/42 + 

+ E 2 (4)f3 2 J- 5 l/2l + J~ 3 1/21 + 2J~ 2 5/42 )], 



R* = AQ* + E(4)f]J- 2 l/3 + 2J~ 1 l/3 + J~ 2 5/42, 



S* = 4J{y/dpQ* + v / ^(S(4)/3J~ 2 (l/3) + E 2 (4)f3 2 ,r 3 (1/3) + £(4)/?J~ 3 (5/42)) 
+ ^/9p{E{A)(3J- 2 {\/3) + E(4)^/6pf3J- 3 (5/42) + V^J~ 3 (1/21))} + 
+ V^{( \$\R* + |e 2 2 | J~ 2 (l/7) +p 2 J" 2 (l/7) ) + (0.7976/3)}. 

The constants A, £lp, [3, \f0~p, \f0p~, \(,q\, an d P2 o,re determined by the par- 
ticular choice of inclusion, while E{4) — I6C2/C5 < 0.7619. 

All bounds obtained here are expressed in terms of the Catalan numbers, area 
fractions and geometric parameters that appear in the Poincare inequality and in 
an extension operator inequality. We start by listing these parameters and give the 
background for their description. It is known Necas (1967) that any H 1 (P) function 
4> can be extended into P as an H 1 (Q) function E(<j>) such that E(<fi) = <p for y in 
P and 



\\E(cj>)\\ Hl(P) <AU\\ H1{P) (5.2) 

where A is a nonnegative constant and is independent of (f> depending only on P. For 
general shapes A can be calculated via numerical solution of a suitable eigenvalue 
problem. Constants of this type appear in Bruno (1991) for high contrast expansions 
of the DC fields inside frequency independent dielectric media. The second constant 
is the Poincare constant Dp given by the reciprocal of the first nonzero Neumann 
eigenvalue of P and we have that Up = 1 + Ep. The last two geometric constants 
appearing in the bounds are the volume fractions Op and Op of the regions P and P. 
Using that C m < 4 m (see section [c]), theorem (5.1l shows that ^2p m n m , ^Pmjf 



and ^2 £m 7 ? m ar e convergent for 77 < 1 /4 J, so that one may prove the following 
theorem 

Theorem 5.2. (Solution of the Eigenvalue Problem) Let R = 1/4J, where J is the 
number prescribed by theorem (5.1l. Then X)m=o?m 7 7 rn converges as a series of real 



numbers and X)m=o V J m r / m converges in the iJ^Q) Sobolev norm for f] < R. More- 
over, u 2 — c 2 k 2 J2m=o {.mV™ an d h„ = h J2m,=r> i m 4>m'n m satisfy the eigenvalue 
problem given by the master system \2.Q\ (or by (|1.3| ). 



Article submitted to Royal Society 



Sub-Wavelength Plasmonic Crystals 



13 



(&) The p m , p m and £~j Inequalities — Stability Estimates 

We now derive the inequalities which bound p m , p m and in terms of lower 
index terms. These inequalities follow from stability estimates for (3.4 3.5 3.6). 

Theorem 5.3. Let m > be an integer. Then 

Pm < ttp[A{2q' m _ 2 + 2q } + 2q' m _ 2 + 2q 

m—o 1 y m—4 



+9 m-4 



+ p m _ 2 + 2 A{2g ;„_ 3 + 2q ' m _ 4 + q ' m _ 5 + q 



+2q' m _ 3 + 2q 

m—4 + 1 m-5 + Q m-5 + Pm-3 

+ 2p 



< -4pm+<7m-2 + 2 P m— 1 ~r Pm—2 

\C-i\ < (M Q 1 {Vo P q / ; n -i + V^pPm + 

+ \/0p(q m-2 + 9 m-3 + 9 m-3) + 

where the p rn inequality holds for m > 2 only. 
Here we have introduced the notation 



(5.3) 



'(9 m-2 + 9 m-3 + 9 m-3)}) 



\8\Pm-l, <I m = Pm-t\$-Mj\ 



1 m-1 



It2|_(«<m-1) 



Pm-^-,11^ 



Proof. We start by proving the p TO inequality Recalling that (3.51 is satisfied by 
V> m in P gives 

A-0m = "0m + G m , in P 
^mlp = 1pm\p, On dP 



where G m = (—i) e t$ip m -2-t - 2k-Vip m _ 1 
sition ip m = u m + v m , where 



ipm-2- Write the orthogonal decompo- 



and 



Au rn = u 

u m = ip. 



Av m = v r , 
v m = 0, 



m , in P 
, on dP 

G m , in P 

on dP. 



(5.4) 



We then have by the triangle inequality that 

IIVwHiiiCP) < ||w to ||hi(P) + ||«m||Hi(P)- 

The term ||wm||j? 1 (p) is bounded using (5.2) and 



to obtain 



\u m \\Hi(P) < \\E(ip m )\\m(p) 



\\u m \\ Hl(P) < A\\lp m \\ Hl{ py 



(5.5) 



(5.6) 



(5.7) 



Here (5.6 1 follows from the fact that the solution of (5.4) minimizes the H 1 (P) 
norm over all functions with the same trace on dP. The term H^mllflifp) can be 
bounded using a direct integration by parts on the BVP for v m 



||«m||ffi(P) < l|Gm|U 2 (P)- 



(5.8) 
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Now, 



\G 



m|| £ 2 (.P) = IK - *) it^m-l-l - 2/t-VVm-l ~ i>m-2\\ L 2 (P) 

< l^ 2 |||V'm-2-£||L 2 (P) + 2|K||||VV' m -2| 2 ||L 2 (P) + 1 1 V4>-2 1 1 L 2 (P) 
m—1 ~r Pm— 2? 



where p m — \\ipm\\Hi(P)- Using (5.7 1 and (5.8 1 in (5.5 1 gives 

Prn _ Ap m + q' rn -2+ 2 P m—1 ~r Prn— 2 j 



(5.9) 



and the p m inequality is established. We now prove the p m inequality. In the weak 
form (3.8 1, set v = ipm in P and v — u m in P to obtain 



+ c r m _ 4 ]u„ l )p + 



+ ([Vf/>m + 



[/c-Vl/>m-l + ^m-2]^m>P = 0, 



We now use the Cauchy-Schwarz inequality on the product of integrals appearing 
in each individual term. For the convolutions, we obtain 



■Vu 



ro/P 



|(V((- l )Vm-2-^£ 2 )-V U „ l ) P | 
l(-*)^£ 2 (VV'm-2-rVw m )p| 



in ' 



where we used that Hit 



m||ffi(P) 



< Ap m . For the double-convolutions, we obtain 



\{v'm-A u m)p\ < ?m-3^Pm- 

Proceeding similarly with the other terms, we obtain 

(VVwW> m ) P < p m (A{2q' m _ 2 + 2q' m _ 3 + q' m _ 4 + q' r ' n _ 4 } + 

m— 3 ~r" 9 m— 4 ~r ^ m— 4 

+ Pm-2()5.10) 

Since the functions "0m have zero average in P, we have the Poincare inequality 



(V4)p < c p (v^ m -vv m )p, 



(5.11) 



where the constant Dp can be computed from the Rayleigh quotient characteriza- 
tion of the first positive eigenvalue for the free membrane problem in P. A simple 



computation using (5.11 1 then gives p^ < Clp{'Vip m -'Viprn)pi where fip = D 2 p + 1. 



Using this inequality (5.10) gives: 

Pm < np(A{2q' m _ 2 + 2q' m _ 3 + q' m _ 4 + q , ; i _ 4 } + 

+^'m-2 + 2<l'm-3 + Q'm-i + q' m -4+Pm-2 + 2p m -l). (5.12) 

It will turn out to be to our advantage to apply (5.121 to the last term 2p TO „ 1 in 



(5.12) so as to replace it with 

1 < n P (A{2q , m _ 3 + 2q> m _ 4 + q , m _ 5 + q'; n _ 5 } + 

1 H ill — o 

+ 2q 

m— 4 Q m— 5 + Qm-5 + Pm-3 

+ 2p m _ 2 ). (5.13) 



Pn 
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Using (5.13) in (5.121 yields the p m inequality in (5.3 > , valid for m > 2 (for to = 1, 
use ( |5.12| ): 

Last we establish the inequality. Setting v = 1 in the weak form (3.- 

obtain 



(/t-Vo-' q - al — o 



a 'm-A/P ' 



(5.14) 



(recall that for to odd, each term on the left-hand side of the above equation vanishes 
individually). Solving for £ 2 „_ 2 we then obtain 



5(-i) m " 2 l 2 (WQ = (^V^- 3 



rn — 3 



°m-4/P 



1-4/ P 



(5.15) 



Pi 



where <j[*_ 2 + (—i) m 2 £ 2 „_ 2 = °m-2- We shall be using this equality for to > 5 
only, so that (c' n *_ 2 )p = and {i/) m -2)p = 0. Moreover, using that (\ip m \)p < 
V^p(IV'«i| 2 )p and (IV'mDp < V^p(|V'm| 2 )p J where 0p and 0p denote the volume 
fractions of the regions P and P, we have that (ip m )p < yf&PPm and (ip m )p < 
\j9ppm- Thus, proceeding with ( 5.15[ ) as we did in the previous stability estimates, 
we obtain 



l&-al < Wo)Q 1 {>/ep~P- 



'm-2-e\tl 



2\l<m-2 



+VdpPm-l+V6p{q' m -3+<lm-4+<l'm-A)- 



Since the iteration scheme at each step involves p m and p m and Cm-l we adjust 
subscripts in the above inequality to obtain the £^-1 inequality 



Id-ll < <^0) Q 1 {y^'Z™-l + V / ^?'m+- 
+ V / «p(?m-2 + C3+4-3)}' 



'(9 m- 2 +Qm-3 + Qms) + 

(5.16) 
□ 



(c) TTie Catalan Numbers 

We briefly present some facts about the Catalan numbers which will be used in 
the sequel. The Catalan numbers C m are defined algebraically through the recursion 



C m +1 — C m -iCl, Cq — 1. 



(5.17) 



These numbers arise in many combinatorial contexts Lando (2003) as well as in 
the study of fluctuations in coin tossing and random walks Feller (1968). It can 
be shown that C m = i 2 ™) and a simple computation then gives the ratio of 
successive Catalan numbers 



a 



m+l 
Cm, 



6 



and 



a 



m+l 



1 

4 



8m + 4 



(5.18) 
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k 





1 


2 


3 


4 


pI 


1 


1/3 


5/42 


1/21 


1/42 



Table 4. Values of p\, where = C m -k/C„ 
The first inequality above provides the exponential bound 

C m < 4 m . 



(5.19) 



It will be convenient to introduce the notation p^ n = C m -k/C m . From (5.18), it 
is clear that p 1 ^ is decreasing in both m and k. In section [d] we shall make use of 
Table 4, in which values of p\ are listed. 

(i) The Even Part of the Catalan Convolution 

The fact that £^ dd = needs to be taken into account in order to provide a 
suitable upper estimate on the incomplete convolution term q' m _i appearing in the 
im-i inequality in (5.3 1. Thus we consider the convolution C n -tCi with the odd 



values of the index i omitted and denote it by C n -iCf even \ We then define the 
even part E(n) by 



E(n) 



Cn-lC, 



(5.20) 



The following lemma gives the estimate E(ri) < E(4), n > 4. 

Lemma 5.1. The following two statements are true for all m > 0: (i) E(2m) is a 
decreasing sequence; and (ii) E(2m + 1) = 1/2. Thus, for all m > 4, we have that 
E(m) < max{E (A), 1/2} = E(4). 

Proof. Statement (ii) is actually just an observation, as one can see by writing out 
the sum C n -iCi. Statement (i) can be deduced from the identity C2m-2iC2i— A m C m 
Koshy (2008). Indeed, dividing both sides of this identity by C2 TO +i, we obtain 



E(2m) = 



1 4CU 4C, 



m+l 



AC 



2m 



4 CVn+l C, 



m+2 



c 



2m+l 



From 5.18 each of the above fractions AC m /C m +e, ^ = 1, 2,..., m+l, is a decreasing 
sequence in m so that their product is also decreasing in m. This completes the 
proof. □ 

(d) Proof of the Catalan Bound 



Proof. (Catalan Bound, theorem (5.1 1) Fix the values of the geometric parameters 
A, Op, Bp and 8p in (5.3 1. Starting with the initial estimates po — OpiPo < 

= 1,2,3,4 



/Wp and |^q| < 1, the inequalities (5.3 1 can be used recursively for 



to determine a number Ji such that 

\e m \<PC m J?, 0<m<4. 
We now proceed inductively: assume that 

Pn i Pn i 

\C\</3C n J n , n G {0,1,2,..., m- 



1}, 



(5.21) 
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where m > 5. We then get for the single convolutions 

q'm-k — Pm-k-e\£i\ 



l-k-l 



)((3CiJ e ) {e even) 



o2 jm—k/~i 

P J ^Tn—k—V^i 







< E{A)(3 2 J 
= E(4)(3J- k 



2 rm—k/^i S~1 



-1-k 



= £J(4)/3J- fe Pm 



where = C m _^/C m and lemma (5.1l was used to introduce the factor £?(4). 
Similarly, for double convolutions we get 



_ fe < E' z (4)f3 2 J- k p^ 2 pJ m G 



where the factor E(4) 2 comes from using lemma (5.1l twice. For the non-convolution 
terms we get p m -k < J k P% l PJ m C m . The same bounds hold for the terms p m -k, 
Qm-k anc ^ Qm-k' so that we have 



Pm — ki Pm — k 

9 m—k 5 y m—k 
- II II 
Q m—k ' 9 m — k 



< 
< 
< 



E{A)f3J~ k p h ~ x $J m C v 



(5.22) 



The proof now essentially consists of applying these bounds to all terms in 
inequalities (5.3 1. The factor J~ k appearing on the right-hand side of each inequality 



is the workhorse of the proof: by taking J sufficiently large, it will allow us to close 
the induction argument. The incomplete convolution term q' m _i presents special 
difficulties, since attempting a bound of the type (5.221 for this term does not 
produce a factor of J~ fc (actually, it produces J° = 1). 
Recall the p m inequality from (5.3 1 



pm < ttp[A{2q , m _ 2 + 2q , m _ 3 + q' m _ A + q" n _ 4 } + 

+ 2 <Zm-2 + 2 9m-3 + <7m-4 + 9m-4 + Pm-2 + 
+ 2 n P (A{2q ;„_ 3 + 2q ' m _ 4 + q ' m _ 5 + q m _ 5 } + 

+ 2 ?' m -3 + 2 9ra-4 + 9m-5 + Qm-5 + Pm-3 + 2p m - 2 )] 



Using (5.221 on this inequality gives 

Pm — Qm @J Cfl 

where Q m is the following polynomial in J -1 



(5.23) 



Q m = Q P [ A{2E(A)(3J- 2 pl + E(4)pJ- 3 p 2 m + E{4)0J-*f^ + E 2 (A)(3 2 J- 4 p 2 m }+ 
+ 2E(A)(3J~ 2 p] n + 2E(A)pj- 3 p 2 m + E(4)pj-*p 3 m + E 2 (A)f3 2 J~ 4 p 2 n + .r 2 p 2 m + 
+ 2n P (A{2E(4)pj- 3 p 2 m + 2E(4)f3J-y m + E(4)f3J- 5 pt + E 2 (4)(3 2 J- 5 p 3 n }+ 
+2E(A)(3J- 3 pl+2E(A)pj- 4 pl+E(A)pj- 5 pt n +E 2 (A)(3 2 J- 5 pl l +J- 3 pl+2J- 2 p 2 n )}. 
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Since we shall be using this inequality for m > 5 only, table [c] can be used to bound 
the numbers p^, so that we may write Q m < Q*, where 

Q* = ftp[ A{2E(A)f3J- 2 l/S+E(4:)/3J- 3 5/A2+E(A)f3J- 4 l/2l+E 2 (4:)/3 2 J- 4 5/A2}- 
+ 2E(4)f3J- 2 l/3 + 2E(4)/3 J~ 3 5/42 + E(4)/3 J~ 4 1/21 + E 2 {4)0 2 J^ 4 5/42 + J~ 2 5/42 
+ 2fl P (A{2E(4)f3J- 3 5/42 + 2E(4)pj-' l l/2l + E(4)pj- 5 l/42+E 2 (4)(3 2 J~ 5 l/2l}+ 
+ 2£(4)/?J- 3 5/42 + 2E{A)(3J- i l/2l + E '(4) /3J- 5 1/42+ 

+ E 2 (4) /3 2 J" 5 1/21 + J~ 3 1/21 + 2J~ 2 5/42 )]. (5.24) 

The strategy now is to determine similar polynomials R m and S m —i for the other 
two inequalities, that is p m < R m (3J m C m and |£ 2 n _i| — S m -i /3</ m_1 C m -i, and 
then take J large enough that all three polynomials are less than unity, allowing 
us to complete the induction argument. Having obtained Q m , it is straightforward 
to obtain R m . Indeed, using (5.22) and (5.231, the p m inequality in (5.3) yields 
Pm < R m /3J m C m , where 

Rm = AQ m + E(4)f3J- 2 p m + 2J- l Pm + J- 2 p 2 m . 

Thus, R m < R* , where 

R* = AQ* + E(A)(3J- 2 l/3 + 2J~ 1 l/3 + J~ 2 5/42. (5.25) 

The £ m _ 1 inequality requires a little more care due to the presence of the incomplete 
convolution term q'£ l _ l . For the remaining terms, we proceed as we did with the 
previous inequalities: 



p\e m - 3 \) < 



VQpPm + + 9m-3 + 9m- 3 ) + V»p(9"L 2 + V«p|^fm-3-<l + 

{VVpQm + ^/8p(0.5pj- 2 pl n + 0.25f3 2 J- 3 pl + 0.5/3J~ 3 p 2 „) 
+ ^(Q.5f3J- 2 pl + 0.5^(3 J- 3 p 2 n + J- 3 p 3 m )}pj m C m . 

since this is an upper bound on we must replace the term (3J m C m with 

PJ^Cm^x as follows: 



(3J m C m = J 



Cm-l 



0J 



m— 1 r 1 



< 4JpJ m ~ x C m - X . 

Using this replacement and the bounds [c] on the numbers pj^, we obtain the upper 
bound 

4J{^9pQ* + \/(h(E{4)f3J- 2 (l/3) + E 2 (4)p 2 ,r 3 (l/3) + E{4) f3 J~ 3 {h / 42)) 
+ V^(#(4)/3J~ 2 (l/3) +£(4) v /^/3J~ 3 (5/42) + ^/OpJ' 3 (1/21))} f3J m - l C m J$.2&) 



It remains to deal with q m _^ To do this, we first write q' m _ x as follows 

Qm-i =p m -i\e \ + Pm - 3 \e 2 \ + P 2\e n -z\ + Pm -i-t\et\ 2<e<m - 3 - 



(5.27) 
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The non-convolution terms then give 



Pm-llfol +Pm-z\&\ +Vl\i 



m — 3 I 



< |^|i? m _ 1 + |^ 2 |J 



_2 Cm— 3 



— 1 C rn _ 1 ^ 

<(|Co 2 |^ + l^ 2 k- 2 (l/7)+p 2 J- 2 (l/7))/3J m - 1 C m _ 1 , (5.28) 



since C m _3/C m _i < C2/C4 = 1/7, if m > 5. The remaining term p m ~i~e\^j\ 2<e<m 3 
is treated in a completely different manner: 



Pm-i-t\t 



2\2<e<m-3 



= Pm~h\£,l\+p m -l\tl \ + ■ ■ ■ +P6\Ci-7\+Pi\^m-5\ 

< {C m -§C± + C m -7Ce + ■ ■ ■ + CeCm-7 + 1 J m 1 

(s~i r<{£ even ) nr< n nr< n \ «2 jm—l 

— y-Jm—l—l^i — ■ il ^2<^m-3 — ^O^m-l) P J 

— {E(4:)C m -i-gCt — 2C 2 C m -3 — 2C C TO _i) 1 J m 1 
= {E{A)C m — 2C2C m _3 — 2CoC TO _i) /3 2 J m 1 

p E(4)C rn - 2C 2 C m - 3 - 2C C m -i \ njm-\Q i 



< /3( J B(4)4-l/4-2)/3J m - 1 C ro _ 1 , 

< (0.7976/3) Pjn-^Crn-l, 



(5.29) 



Thus, adding (5.26), (5.281 and (5.29), we set 



S* = AJ{^f¥pQ* + V^(^(4)/3J" 2 (l/3) + E 2 {A){3 2 J~ 3 (l/3) + £(4)/3J~ 3 (5/42)) 
+ V^(£:(4)/3J- 2 (l/3) + £(4) 7<9p/3J- 3 (5/42) + y/ffaj-* {1/21))} + 

+ V^{( \$\R* + \&\J- 2 {l/f) + P2 J- 2 (l/7) ) + (0.7976/3)X5.30) 

Thus, taking J 2 such that Q*(J 2 ), R*(J 2 ), S*(J 2 ) < 1 and J = max{Ji, J 2 }, we 
have shown that the induction hypothesis (5.211 implies 



Pm, Pm, \£,m\ ^ f3C m J m , 



so that in fact (5.31 1 holds for every integer m. 



(5.31) 
□ 



(e) Proof of Theorem 5.2' Solution of the Eigenvalue Problem 



Proof. The weak form of the master system is 

e-\V + i V k)h v (y) ■ (V - ir)k)v(y) - if %h^)v{y)\ = for all v e H^JQ). 

(5.32) 

Using that e" 1 = lf&/{rf'^. — 1) in P and e v = 1 in P, and multiplying by 
(v 2 £, 2 ~ 1); gives the equivalent system 



a v (h, C 2 ; v) = for all « G H^XQ), 
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in which 



a v (h, £ 2 ; v) = — j (V + irjk)h ■ (V — irjk)v + 

+ / [v 2 £ 2 (V + ir]k)h- (\7 -i v k)v- (r] 2 £, 2 -l)ri 2 £ 2 hv] . 
JQ 

This form can be expanded in powers of rj, 

a n {h, £ 2 ; v) = a (h; v) - irja^h; v) - rj 2 a 2 (h, £, 2 ;v) + irj 3 a 3 (h, £ 2 ; v) + 7] 4 a 4 (h, £, 2 ;v) 

(5.33) 

in which the 

Q"fn real forms 

a (h;v) = - j Vh- Vv, 
a\(h;v) = J (hk • VC — V/i • kv), 
a 2 (h,t 2 ;v)= [ hv-Z 2 [ (Vh-Vv + hv), 

JP JQ 

a 3 (h,£_ 2 ;v) = £ 2 / (hk ■ Vv - V/i • kv), 

JQ 

a 4 (h,e;v) = e f (i-e)hv. 

Jq 

Define the partial sums 



e,' N =J2v m e m and h»=f2v m h m . 

m— m— o 

For rj < R, the sequence {£, 2 ' N } converges to a number £ 2 and the sequence {h™} 
converges in Hi to a function h v ; thus 

w) - for all v e ^(Q) and i = 0, . . . , 4. 

Therefore, a,j(h v , £ 2 ; t>), j = 1, . . . , 4, has a convergent series representation in pow- 
ers of rj, in which the m th coefficient is related to the coefficients & and /i^ by 

(j = °) J - Vh ™ ■ 

(j = 1) ^ (ft m K • Vv - Vft m • kv) , 

(.] = 2) / h m v - I (V(C £ 2 /i m -£) • VC + ($h m - e )v) , 

JP JQ 
0'=3) / (($h m -e)k • VC - V(^ 2 /l m _ £ ) • KU) , 

(i = 4) / fe 2 ^ - ijihjh^) v. 

JQ 
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From these, one obtains the m 



coefficient of a TI (h v ,^;v) (see 5.331, which, by 
7 the relations h m = h i m ip m , ($h m -i = h i m a' m and ^j_jh m ^e = 
is seen to be equal to the —i m ho times the right-hand side of equa- 
tion (|3.8|. All these coefficients are therefore equal to zero, and we conclude that 

together with the frequency 

□ 



a ri(h v , v) = 0. This proves that the function h 



solve the weak form (5.321 of the master system. 



6. Effective Properties, Error Bounds and the Dispersion 

Relation 

In this section we start by identifying an effective property directly from the dis- 
persion relation. We then discuss the relation between effective properties and qua- 
sistatic properties. Next we provide explicit error bounds for finite-term approxima- 
tions to the first branch of the dispersion relation for nonzero values of r\. The error 
bounds show that numerical computation of the first two terms of the power series 
delivers an accurate and inexpensive numerical method for calculating dispersion 
relations for sub-wavelength plasmonic crystals. 



(a) The Effective Index of Refraction - Quasistatic Properties and 

Homogenization 

The identification of an effective index of refraction valid for r\ > follows 
directly from the dispersion relation given by the series for ljz. Indeed the effective 
refractive index n^ s is defined by expressing the dispersion relation as 

< = ^ (6-D 

"eff 

and it then follows from the expansion for u>^ that the effective refractive index has 
the convergent power series expansion 

oo 

„-2 _ „-2 , „2m f 2 



= n 



<|S 



(6-2) 



We now discuss the relationship between the effective index of refraction and the 
quasistatic effective properties seen in the d — > limit with k fixed. The effective 
refraction index n c s can be rewritten in the equivalent form by the equation n^ s — 



By setting v = h v in the weak form of the master system (5.32), it is easily 
seen that ^ > for all rj within the radius of convergence, so that n^ ff > for 
those values of r\. Following Pendry et al. (1999), see also Kohn & Shipman (2008), 
we define the effective permeability by 

(B 3 ) cS 

Meff = fTT s » (6-3) 

and we then define e e g through the equation 

n eff = EeffMefif- (6.4) 
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The quasi-static effective properties are recovered by passing to the limits 



,2 u„ „2 

qs 



= lim n cS , /j qs = lim /j off , e qs = lim e cS . 

17— >0 77^0 r\— >0 



A simple computation shows that fj, qs = (V'o)q > ( see appendix). Hence, we have 
that )i c g > for 77 in a neighborhood of the origin, so that e e g > for these values 
of 77, since n 2 s > for all 77 in the radius of convergence. Thus, one has a solid basis 
on which to assert that plasmonic crystals function as materials of positive index 
of refraction in which both the effective permittivity and permeability are positive. 

For circular inclusions we have used the program COMSOL to compute that 
("0o)q ~ 0.98, so that only a mild effective magnetic permeability arises. 

Having established that h v (y)e lr > K ' y is the solution to the unit cell problem, we 
can undo the change of variable y = fcx/77 to see that the function 

h v e i(*fr*-«,*) ) ( 6 .5) 

where h v is the Q-periodic extension of h v to all of M 2 , is a solution of 

V x • (e^Vx/g = ^d tt h v , xel 2 , (6.6) 

for every 77 in the radius of convergence. 

We investegate the quasistatic limit directly using the power series (6.51. Here 
we wish to describe the average field as d — > 0. To do this we introduce the three- 
dimensional period cell for the crystal [0, d] 3 . The base of the cell in the X1X2 plane 
is denoted by Qd — [0,d] 2 and is the period of the crystal in the plane transverse 
to the rods. We apply the definition of B e g and H e g given in Pendry et al. (1999) 
which in our context is 

(B 3 ) eS =^JK (^j e^-^d Xl dx 2 (6.7) 



and 



1 KO,o,d) ^ 

(ff 3 )eff = L J M — ) e i(M - x -"'*>dx3. (6-8) 



^ J (0,0,0) V V J 



Taking limits for k fixed and d — > in (6.5 1 gives 

lim(5 3 )cff - (^Qhoe 1 ^-^ and rim(tf 3 )eS = h^^-^ , 

9 2 k 2 

in which uj„ a — These are the same average fields that would be seen in a 

qs n 2 3 t> 

quasistatic magnetically active effective medium with index of refraction n qs and 
;Uq S that supports the plane waves 

(B 3 ) qs = l i qa h e i ^ x - u ^ and (H 3 ) qs = j^C**-*—,.*), 

where fi qs — (V'o)q. It is evident that these fields are solutions of the homogenized 

n 2 

equation -^dttu = Au. This quasistatic interpretation provides further motivation 
for the definition of n B g for nonzero 77 given by |6.2| 
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Now we apply the definition of effective permeability fj, c g given in Pendry et 
al. (1999), together with n e g to define an effective permeability e c s for 77 > 0. The 
relationships between the effective properties and quasistatic effective properties 
are used to show that plasmonic crystals function as meta-materials of positive 
index of refraction in which both the effective permittivity and permeability are 
positive for r\ > 0. 

(&) Absolute Error Bounds 

The Catalan bound provides simple estimates on the size of the tails for the 
series £ 2 and h v , 



m-mo + 1 



2 and E mo , 



h ^f 

m— mo + 1 



We have established convergence for 77 < 1/4 J, so that we may write 77 = a/4 J, 
< a < 1. Then, using that |f§J < [iC 2m J 2m , C 2m < 4 2m and AJ-q = a, we have 



\E, 



Similarly, for h v we have that 



< 



V F 2 n 2m 



m— mo + 1 

00 



< P J2 C 2mJ 2m ri 2m 



m— mo + 1 
00 



m— mo + 1 
00 

Yl « 2m 



m— mo + 1 
3 ,2m +2 

I -a 2 ' 



\\E mo , h \\m( Q ) < 2f3\h \ Y 



(6.9) 



(6.10) 



(c) Relative Error Bounds 



In this section, we use the absolute error bound (6.10) with mo = 1 to obtain a 



relative error bound for the particular case of a circular inclusion of radius r = 0.45 
Shvets & Urzhumov (2004). The first term approximation to is 



h v = h tp + ihoipiT) + E 1<h . 



(6.11) 



For a circular inclusion of radius r = 0.45, we have J < 85, (3 < 0.79, ||V'o|| = 
0.97 and = 0.02, where || • || = || ■ \\h 1 (q)- Thus, using bound (6.10 1, the relative 
error R± h is bounded by 



\Ri,h\ = 



\Ei, h \ 



1.58, 



1.58 



\h + ihoipiriW 



< 



IM - ||ftoVi|||»7l 



< 



l-Q 



0.97- 0.02^ 
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Figure 3. Graph of the first branch of the 
Figure 2. Solid curve is Ri i6 and dotdash dispersion relation. Dashed curves represent 
curve is Ri,h- error bars. 

so that for a < 0.2 the relative error is less than 8.2%. The graphs of ipo and ty\ 
can be found in the Appendix. The first term approximation to ^ is 



(6.12) 



In the Appendix we indicate how the tensors may be computed. For an inclusion 
of radius r = 0.45, we have £q = 0.36 and £f = —0.14. Thus, using bound ( 6.9 1 , the 
relative error is bounded by 



\Ri„ 



\Ei. 



< 



" 1 — OL 



< 



0.79 



l-Q 2 



l^o + C 2 VI - l£ 2 + £ 2 VI 

so that for a < 0.3 the relative error is less than 2% 



0.36 - 0.143^ 
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Appendix A. Explicit Expressions for Tensors 



The tensors £q and £| were calculated using the weak form (3.8 1, as follows. Setting 
to = 2 and v = 1 in (3.8 1 and solving for £q gives 



2 _ (ft -S/tpi + il) )p 



(V>o)c 



Setting 777 = 4 and v = 1 and solving for gives 



Q 



(Al) 



(A 2) 



All integrals appearing in (Al) and (A2| were then computed using the program 
COMSOL. 



Appendix B. Bounds on p , p and |£q| 

We have that po = Op, po < \jQp, and |£q| < 1- These bounds are obtained as 
follows: we have ipo = 1 in P, so that po = 8p. Using the BVP for ipo in P one can 
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prove that < V'o(y) < 1, Vy G P, and that p\ = (ipo)p. These two facts together 
give the estimate po < y/Op. Setting v = ip\ in the weak form for ipi, we get that 
(k-\/tpi)p = — (V^i ■ Vt/>i)p < 0. Using this in expression (All gives £ 2 < 1. Since 
Op, 6p < 1, these three estimates allow us to take (3 < 1. 



Appendix C. Computing the Constant A for Circular 

Inclusions 

Given a function ip € Hp er (P), let u S H 1 (P) satisfy 

V 2 w - u = in P, 
u = ip on 9P. 

We seek to compute a number A such that ||w||jji(p) < ^IIV'llj^p) f° r an "0- 
Following Bruno (1991), we will calculate a value of A for circular inclusions P of 
radius r < 0.5 by restricting ip to the annulus between P and the circle of unit 
radius. It suffices to consider real- valued functions ip that minimize the H 1 norm in 
the annulus, that is V 2 ip — ip = 0, ro < r < 0.5. A function of this type is given 
generally by the real part of an expansion ip(r, 9) = Y^=o ( c «^n( r ) + d n K n {r)) e vn6 , 
in which c n and d n are complex numbers and /„ and K n are the "modified" Bessel 
functions. The continuous continuation of ip into the disk with V 2 it — u = is given 
by the real part of u(r, 9) — X^Lo fnln(x)e under the relations 

fn=c n + d n §p\. (CI) 
K n {r ) 

One computes that 1 1 R.e "0 1 1 2 = illV'll 2 an d ||Reu|| 2 = |||u|| 2 , so we may work 
with the complex functions rather than their real parts. The Hclmholtz equation 
in P and integration by parts yield 

Ml 2 



HHP)= / Vw | 2 + H 2 )dA= / ud n u, 
Jp JdP 

and this provides the representation IHI^i = 2nri J2^=o fnln( r i)fnl'n( r i)- The 
analogous representation in the annulus is 

\\4>\\h= r(dMl,9)W^jd9- r (dMro,0)W^)d9 
Jo Jo 

oo 

= 2tt J2 (c«/n(0.5) + d n K n (0.5)) (c n l' n (0.5) + d n K' n (0.5)) + 

n=0 

oo 

- 2irr 2J (c„I„(r ) + d n K n (r )) (c n I' n (r ) + d n K' n (r )) 

We seek a positive number A such that, for all choices of complex numbers 
c n and d n we have < ^4||t/>|| 2 — \\u\\ 2 . The right-hand-side of this inequality is a 
quadratic form in all of the coefficients (c n ,d n ) that depends on A, 

A\\ip\\ 2 - \\u\\ 2 = m"c„c„ + m 12 c n d n + vn^d n c n + m 22 d n d n , 
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in which the depend on A and are conveniently expressed in terms of the 
functions 

II n (r) = I n (r)i;(r), KK„(r) = K n (r)K' n (r), IK„(r) = I n {r)K' n {r) 
KI„(r) = tf„(r)/;(r), JJ n (r) = W 

= -r H„(r ) -Ar II„(ro)+AII„(0.5), 

mf = -r JJ„(r )- Ar KK n (r ) + AKK(0.5), 

m 12 = - ro KI(r )- Ar KI n (r ) + AKI„(0.5), 

ml 1 = -r KI(r )- J 4r IK n (r ) + ^IK n (0.5). 

The form is Hermitian, as one can show that m^ 2 = to 21 by using the fact that 
rWron[/ n , K n ] is constant. 

We must find A > such that to" > and m^m 22 - m^ 2 m 21 > for all 
n = 0,1,2,.... These quantities are equal to 

m)} = /3(r ,0.5)A-a(r ), 
m^mf-m^ml 1 = e(r , 0.5)A 2 - <5(r , 0.5) A, 

in which 

a n (r) = rll n (r), 
Pn(r,s) = ll n (s) - rll„(r), 

<5„(r, s) = rs [II n (r)KK„(s) + II„(s)JJ n (r) - KI„(r)IK„(s) - KI„(s)KI„(r)] + 

+ r 2 [-II„(r)KK„(r) - II„(r)JJ(r) + KI n (r)IK n (r) + KI n (r)KI„(r)] , 
e n (r, s) = rs [-II n (r)KK„(s) - II„(s)KK„(r) + KI„(r)IK n (s) + KI„(s)IK„(r)] . 

The numbers a„(r) and /3 ra (V, s) for r < s are positive; the latter because (rlnl^)' = 
-(r 2 + n 2 )l 2 + rl' rl > 0. One can show that 6 and e are positive. Thus, it is sufficient 
to find A > such that, for all n = 0, 1, 2, ... , 



^4 > max 



a n (ro) 6n(Vo,0.5) 



/3 n (r ,0.5) ' e„(r ,0.5) 
Table [Appendix D] shows computed values of A for various values of ro- 

Appendix D. Graphs of ip m and Table of A, ftp and J for 

Circular Inclusions 



Article submitted to Royal Society 



Sub-Wavelength Plasmonic Crystals 



27 




Figure 4. Graph of tpo- This function is Figure 5. Graph of tpi when k — (1,0). This 
symmetric about the origin. function is antisymmetric about the origin. 



r 


0.1 


0.2 


0.3 


0.4 


0.45 


A 


1.058 


1.293 


1.907 


3.956 


4.840 


dp 


1.0 


1.0 


1.0 


1.0 


1.0 


J 


15 


17 


22 


29 


85 



Table 5. Values of A, ftp and J for circular inclusions of radii rd. 

References 

Auriault, J. L. & Ene, H. Macroscopic modelling of heat transfer in composites with 
interfacial thermal barrier. Internal. J. Heat Mass Transfer, 37:2885-2892, 1994. 

Bakhvalov, N. S. & Panasenko G. P. Homogenization: Averaging Processes in Periodic 
Media. Kluwer, Dordrecht/Boston/London, 1989. 

Bensoussan, A., Lions, J. L. & Papanicolaou, G. Asymptotic Analysis for Periodic Struc- 
tures, volume 5 of Studies in Mathematics and its Applications. North-Holland Publish- 
ing Co., 1978. 

Bouchitte, G. & Felbacq, D. Homogenization near resonances and artificial magnetism 

from dielectrics. C. R. Acad. Sci. Pans, 1(339) :377-382, 2004. 
Bouchitte, G. & Felbacq, D. Theory of mesoscopic magnetism in photonic crystals. Phys. 

Rev. Lett, 2005. 

Bouchitte, G. & Schweizer, B. Homogenization of maxwell's equations with split rings. 
Preprint, 2008. 

Bouchitte, G. & Felbacq, D. Left-handed media and homogenization of photonic crystals. 
Optics Letters, 30(10):1189-1191, 2005. 

Bouchitte, G. & Felbacq, D. Negative refraction in periodic and random photonic crystals. 
New Journal of Physics, 159, 2005. 

Bruno, O. The effective conductivity of strongly heterogeneous composites. Proc. Roy. 
Soc. Lond. A , 433(1888):353-381, 1991. 

Cherednichenko, K., Smyshlyaev, V. & Zhikov, V. Non-local homogenized limits for com- 
posite media with highly anisotropic fibres. Proc. R. Soc. Edinburgh, 136A-.87-1 14, 
2006. 

Cioranescu, D. & Paulin, J. Homogenization in open sets with holes. J. Math. Anal. Appl., 
71(2):590-607, 1979. 

Conca, C.,Orive, R. & Vanninathan, M. On burnett coefficients in periodic media. Journal 

of Mathematical Physics, 47(032902):1-11, 2006. 
Donato, S. & Patrizia, M. Monsurro. Homogenization of two heat conductors with an 

interfacial contact resistance. Anal. Appl. (Smgap.), 2:247-273, 2004. 
Feller, W. An Introduction to Probability Theory and its Applications, volume 1. Wiley, 

3rd edition, 1968. 

Jackson, J. Classical Electrodynamics. John Wiley & Sons, Inc., third edition edition, 1999. 



Article submitted to Royal Society 



28 



Fortes, Lipton & Shipman 



Kamotski, V., Matthies, M. & Smyshlyaev, V. Exponential homogenization of linear sec- 
ond order elliptic pdes with periodic coefficients. SIAM Journal on Mathematical Anal- 
ysis, 38(5):1565-1587, 2007. 

Kohn, R. & Shipman, S. Magnetism and homogenization of micro-resonators. SIAM Mul- 
tiscale Model Simul, 7(l):62-92, 2008. 

Koshy, T. Catalan Numbers with Applications. Oxford University Press, 2008. 

Lando, S. Lectures on Generating Functions, Stud. Math. Lib. 23 American Mathematical 
Society, 2003. 

Lipton, R. Heat conduction in fine scale mixtures with interfacial contact resistance. SIAM 
J. Appl. Math, 58(l):55-72, 1998. 

McPhedran, R.C., Poulton, C.G., Nicorovici, N.A. & Movchan, A. Low frequency Cor- 
rections to the Static Effective Dielectric Constant of a Two Dimensional Composite 
Material, Proc. Roy. Soc. A, 452, 2231-2245 (1996). 

Necas, J. Les methodes directes en theorie des equations elliptiques. Masson et Cie, 
Editeurs, Paris, 1967. 

Odeh, F. & Keller, J. Partial differential equations with periodic coefficients and bloch 
waves in crystals. Journal of Mathematical Physics, 5:1499-1504, 1964. 

Panasenko, G. Boundary conditions for the high order homogenized equation: laminated 
rods, plates and composites. C. R. Mecanique, 337:18-24, 2009. 

Pendry, J., Holden, A., Robbins, D. & Stewart, W. Low frequency plasmons in thin-wire 
structures. J. Phys.: Condens. Matter, 10:4785-4809, 1998. 

Pendry, J., Holden, A., Robbins, D. & Stewart, W. Magnetism from conductors and en- 
hanced nonlinear phenomena. IEEE Trans. Microw. Theory Tech., 47(ll):2075-2084, 
1999. 

Povinelli, M., Johnson, S., Joannopoulos, J. & Pendry, J. Toward photonic-crystal meta- 
materials: Creating magnetic emitters in photonic crystals. Applied Physics Letters, 
82:1069-1071, 2009. 

Shvets, G. & Urzhumov, Y. Engineering the electromagnetic properties of periodic nanos- 
tructures using electrostatic resonances. Phys. Rev. Lett., 93(24):243902-l-4, December 
2004. 

Shvets, G. & Urzhumov, Y. Electric and magnetic properties of sub-wavelength plasmonic 

crystals. J. Opt. A: Pure Appl. Opt., 7:S23-S31, 2005. 
Smith, D., Padilla, W., Vier, D., Nemat-Nasser, S. & Schultz, S. Composite medium with 

simultaneously negative permeability and permittivity. Phys. Rev. Lett., 84(18):4184- 

4187, May 2000. 

Smyshlyaev, V. & Cherednichenko, K.. On rigorous derivation of strain gradient effects in 
the overall behaviour of periodic heterogeneous media. J. Mech. Phys. Solids, 48:1325- 
1357, 2000. 

Smyshlyaev, V. Propagation and localization of elastic waves in highly anisotropic periodic 
composites via two-scale homogenization. Mech. Mater., 41:434-447, 2009. 

Figure Captions. 

1 Unit cell with plasmonic inclusion. 

2 Solid curve is R\£ and dotdash curve is R\,h- 

3 Graph of tp . This function is symmetric about the origin. 

4 Graph of ip l when k = (1, 0). This function is antisymmetric about the origin. 
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